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INTRODUCTION 


The  relationship  between  free-air  gravity  anomalies  and 
sea-floor  bathymetry  is  one  of  the  most  important  problems  in 
marine  gravity.  For  purely  practical  reasons,  deciphering 
the  nature  of  such  a relationship  is  vital  to  predicting 
quantitatively  the  gravity  field  in  oceanic  areas  in  which 
only  bathvmetry  is  known.  Closely  related  to  this  practical 
problem  is  the  question  of  the  physical  basis  for  any  observed 
or  proposed  dependence  of  gravity  upon  topography.  Such  a 
question  bears  on  the  mechanisms  for  creation  and  evolution 
of  oceanic  crust  and  lithosphere  and  on  the  possible  interaction 
of  the  lithosphere  with  motions  in  the  asthenosphere. 

Marine  gravity  and  sea-floor  topography  cannot  be  related 
by  a simple  mathematical  expression  (e.g.,  a line)  that  has 
validity  in  all  oceanic  environments.  A recent  synthesis 
of  submarine  gravity  data  and  bathymetry  by  Woolard  and 
Daugherty  (1970)  demonstrated  the  necessity  to  divide  oceanic 
regions  by  tectonic  types  and  the  difficulty  in  simply 
relating  free-air  anomaly  and  depth  even  within  groups  of 
tectonically  similar  environments.  Degree  averages  of  marine 
gravity  anomalies  and  depths  within  selected  regions  were  also 
correlated  with  varying  degrees  of  success  by  Watts  and  Talwani 
(1974)  , Sclater  et_  al^.  (1975)  , and  Watts  (1976)  . 

A necessary  step  in  improving  what  heretofore  have 
generally  been  strictly  empirical  attempts  to  derive  a rule 
relating  gravity  and  bathymetry  over  a selected  region  is  to 
take  into  proper  account  the  findings  of  plate  tectonic  studies 
of  the  oceans.  At  mid-ocean  ridges,  the  lithosphere  is  thin, 
perhaps  no  thicker  than  the  crust  (Francis  and  Porter,  1973; 
Solomon  and  Julian,  1974;  Orcutt  et_  al . , 1975;  Rosendahl  et  al . , 

1976)  , and  sea  floor  topography  except  at  the  shortest  wave- 
lengths is  isostatically  compensated  with  a shallow  compensation 
depth  (Dorman,  1975;  McKenzie  and  Bowin,  1976;  Watts  and  Cochran, 

1977)  . Sea  floor  topography  created  by  mid-plate  volcanic 
activity  is  compensated  much  deeper  because  of  the  substantially 


3 


thicker  lithosphere,  and  the  compensation  is  regional  (Veniny 
Meinesz,  1941) , involving  flexure  of  the  oceanic  plate 
(Walcott,  1970;  Watts  and  Cochran,  1974;  Watts  et  al . , 1975), 
the  effective  elastic  thickness  of  which  increases  with  sea 
floor  age  (Caldwell  and  Turcotte,  1978) . At  trench  systems, 
topography  is  not  isostatically  compensated  owing  to  the 
large  dynamic  forces  associated  with  lithosphere  subduction 
(e.g.,  Vening  Meinesz,  1954).  Still  controversial  are 
explanations  for  frequent  correlations  in  very  long  wave- 
length gravity  and  topographic  anomalies  in  the  oceans 
(Menard,  1973;  Anderson  et  al^.  , 1973;  Weissel  and  Kayes, 

1974;  Sclater  et  a_l . , 1975;  Marsh  and  Marsh,  1976;  Cochran 
and  Talwani,  1977;  Watts,  1978),  with  asthenospheric  flow 
the  most  exciting  but  still  unproven  hypothesis. 

This  report  lays  the  mathematical  framework  for  a 
simple  but  general  relation  between  topography  and  gravity 
in  two  dimensions  for  stable  ocean  basins.  The  conceptual 
framework  for  the  relation  is  based  on  flexure  theory  for 
thin  elastic  plates  loaded  from  above.  The  evolution  with 
sea-floor  age  of  lithospheric  temperature  and  rheology 
is  abstracted  to  an  effective  elastic  layer  thickness  which 
grows  with  plate  age.  In  order  to  derive  gravity  from 
topography  using  this  relation,  both  the  lithospheric  age 
and  the  age  of  any  more  recently  superposed  volcanic 
constructs  (islands,  seamounts,  aseismic  ridges)  must  be 
known.  Guidelines  are  given  for  estimating  these  ages. 

Finally,  the  quantification  of  the  bathymetry-gravity 
relation  using  data  from  the  central  Pacific  basin,  work 
being  conducted  during  the  second  year  of  this  contract  research, 
is  briefly  described. 
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THEORY 


We  seek  a relation  between  bathymetry  and  gravity  in 
stable  ocean  basins.  We  shall  work  explicitly  in  spatial 
dimensions,  rather  than  use  one-  or  two-dimensional 
wavenumber  representations,  so  as  to  permit  consideration 
of  regions  of  arbitrary  geometry  and  spatial  extent. 

We  regard  the  sea- floor  topography  as  composed  of 
three  parts:  (1)  the  long-wavelength  deepening  of  the 

sea  floor  with  age  due  to  thermal  contraction  of  the 
lithosphere  (Sclater  et  al.,  1971);  (2)  the  volcanic 

topography  emplaced  as  a load  on  top  of  the  lithosphere; 
and  (3)  the  lithospheric  response  to  that  load.  We  assume 
that  the  effect  of  (1)  on  the  bathymetry  may  be  removed 
with  a suitable  age-depth  relation  (Sclater  et  al.,  1975; 
Cochran  and  Talwani,  1977);  its  effect  on  gravity  is 
negligible  far  from  ridges  (Lambeck,  1972). 

Let  the  lithosphere  be  modeled  as  a thin,  spherical, 

elastic  shell  of  thickness  T,  Young's  modulus  E,  and 

Poisson's  ratio  v,  overlying  a fluid  interior  of  density 

p . Let  R be  the  radius  to  the  midplane  of  the  shell, 
m r 

let  g be  the  gravitational  acceleration  at  radius  r = R - T/2, 
and  let  D be  the  flexural  rigidity  of  the  shell 


n = ET3 
u 12(1  - v2) 


(1) 


Consider  a uniform  vertical  load  q (force  per  unit  area) 
applied  to  a circular  unit  area  on  the  surface  of  the  shell. 

The  vertical  deflection  (negative  if  downward)  of  the  lithosphere 
is  given  by 

w = kei  £ (2) 

2 it  (ET/R2  + Ap  g) 

(Brotchie  and  Sylvester,  1969;  Brotchie,  1971),  where  £ is  the 
distance  from  the  load  center,  normalized  by  the  radius  of 
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relative  stiffness 


l 


| ET/R2 


D 

+ 


Ap  cj 


(3) 


kei  is  a Bessel-Kelvin  function  of  order  zero  (Abramowitz 
and  Stegun,  1964) ; and  Ap=  p - p where  p is  the  density 
of  seawater  if  the  load  is  applied  at  the  sea  bottom  and 
zero  if  the  load  is  applied  on  land.  Note  that  (2)  is 
equivalent  to  the  expression  (Brotchie  and  Sylvester,  3969} 
for  the  response  to  a point  force  p since  p = qnd'fc2  where 
^ = d is  the  radius  of  the  circular  unit  area.  Equations 
(2)  and  (3)  may  be  simplified  because  for  oceanic 
lithosphere  ET/R2  ~ 101  and  Ap  g - 2 x 10 3 in  c.g.s.  units, 
so  Ap  q >>  ET/R2  and 


w = 


q 

2itAp  g 


kei  £. 


(4) 


We  may  generalize  (4)  to  a distributed  load.  Let 
q(x,y)  be  the  force  per  unit  area  exerted  on  the  lithosphere 
by  topography.  Then  the  deflection  w is  given  by 


w ( x , y ) 


1 

2irAp  g 


q(x'  , y ' ) 


kei 


( r'/i ) dx'dy' 


(5) 


i , 

where  r'  = [ ( x ’ - x) 2 + (y'  - y)2]  ' 2 . Equation  (5)  may 

be  thought  of  as  the  double  convolution  of  q with  the  flexural 
response  function 


$ (x,y) 


1 

2nAp  g 


kei 


(x2 


(6) 


which  in  turn  is  a function  of  the  local  flexural  ridigity  D, 
or  equivalently  the  effective  elastic  lithospheric  thickness 
T,  through 
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A single  convolution  of  (two-dimensional)  load  with  response 
function  was  also  used  by  Roufosse  and  Parsons  (1977)  in  their 
study  of  the  Hawaiian  ridge.  We  regard  £ as  a function  of 
x and  y,  dependent  in  ? predictable  fashion  on  the  age  of  the 
lithosphere  (Caldwell  and  Turcotte,  1978)  and  on  the  age 
of  the  lo'-al  load.  If  the  load  was  emplaced  in  a time  short 
compared  to  the  lithosphere  age  and  if  viscous  relaxation 
of  stress  subsequent  to  initial  lithospheric  flexure  has 
been  minimal,  then  l is  simply  a function  of  the  time 
difference  between  plate  age  and  load  emplacement  age. 

Unfortunately,  the  load  q(x,y)  is  not  known  a priori 
but  rather  the  ocean  floor  elevation  h(x,y)  with  respect  to 
some  arbitrary  datum  is  known.  As  noted  above,  h contains 
contributions  from  both  the  topographic  load  and  the  litho- 
spheric response.  We  may  resolve  this  difficulty,  however, 
by  iteration.  Assume  initially  that  q(x,y)  = Ap’gh(x,y) 
where  Ap1  = p - and  p^  is  the  density  of  the  major  units 
constituting  the  topographic  variations.  Then  use  (5)  to 
calculate  w(x,y).  Now  set  q(x,y)  = Ap’g(h  - w) . Recalculate 

w from  (5),  and  repeat  the  last  two  steps  until  q and  w 

converge  to  a steady  solution. 

Once  a self-consistent  decomposition  of  h into  load 
topography  q/Ap'g  plus  plate  deflection  w is  achieved, 
the  calculation  of  gravity  is  straightforward.  Two  terms 
contribute  to  the  gravity:  (i)  the  attraction  of  the  topog- 
raphy and  (ii)  the  deflection  of  the  Moho  and  of  any  other 
density  contrast  interfaces  within  the  lithosphere.  Any 
density  contrast  between  the  base  of  the  lithosphere  and  the 
asthenosphere  would  also  contribute  to  (ii) , but  this  contri- 
bution is  probably  negligible.  For  both  (i)  and  (ii)  the 

gravity  anomaly  may  be  written  in  the  form  of  that  due  to 
a surface  mass  distribution  o(x,y) 
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Ag(x,y,z) 


Gz 


dx 'dy ' (7) 


(f 


[ (x- 


O ( X 1 , y 1 ) 

- x) z + Ty1  - y)  r~+  ztT5/2 


where  G is  the  gravitational  constant  and  z>  0 is  the  vertical 
distance  between  the  observation  point  (x,y,z)  and  the 
horizontal  plane  on  which  a is  evaluated.  For  the  topographic 
contribution,  o = Ap'  h on  the  plane  of  the  topographic  datum 
(h  = 0) . For  the  plate  deflection  contribution,  a - Ap  w on 
the  plane  corresponding  to  mean  Moho  depth. 

An  independent  check  on  any  bathymetry-gravity  model  is 
the  lithospheric  stress.  The  stress  field,  as  inferred 
from  lithospheric  earthquake  mechanisms,  say,  may  be  a strong 
function  of  the  load  (Rogers  and  Endo,  1977) . The  stress 
may  also  be  calculated  directly  once  the  load  q(x,y)  is 
determined.  For  the  single  load  q applied  to  a circular 
unit  area,  the  bending  moments  M and  stress  resultants 
N expressed  as  components  along  the  polar  angle  <J>  and 
azimuthal  angle  0 with  respect  to  the  load  center  are 
(Brotchie,  1971) 


M<(>  = ~ 2tt’  (ker  + ^T~  kei'O 


M„ 


= - 57  (^ker  C + kei’C) 


N , = 


N„  = 


qft  2ET/R 


(8) 


(-  + ker’  5 ) 


2 TT  D 

_ ETw  _ - qET/R  [ l2  .1 

4>  “ R 2 rr  I CD  ^ 


+ ker'  5)  + kei  £ 


] 


where  ker  is  another  Bessel-Kelvin  function  of  order  zero  and 
where  prime  denotes  first  derivative. 
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The  horizontal  stress  components  within  the  plate  at  radius 
r are  given  by 


0 v = 


N , 

_ <P 


12M. 

<J> 

~ srr-  (r 


^ a rn 


0 + 12M0 
T U 


- R) 


R) 


(9) 


For  a distributed  load,  q(x, y)  , suppose  0 is  measured  clock- 
wise from  the  +y  direction,  and  let  the  x and  y axes  point 
E and  N,  respectively.  Then 


r 

“ 

0 

XX 

a 

xy 

(x,y)  = 

f 

R 

O. 

<p 

0 

a 

yx 

a 

yy  _ 

— 

0 

a0 

where 


R = 


0 = tan 


sin0  cos0 
cos0  - sin0 

-1  x - x' 


(x'  - x,  y'  - y)  R dx'dy 


(10) 


y - y' 

[ (x1  - x) 2 + (y'  - y ) 2 1 ^ 


R 


and  where  and  Og  are  determined  from  q using  (8)  and  (9) 
with  £ = R6/f. 


SOME  PRACTICAL  CONSIDERATIONS 


The  theoretical  relationship  between  topography  and 
gravity  outlined  above  can  be  applied  to  any  finite  area  over 
which  bathymetry  and  age  are  specified  as  functions  of  position. 
We  are  currently  applying  these  concepts  to  regions  in  which 
degree  averages  of  ocean  floor  depth  (corrected  for  mean 
age)  and  free  air  gravity  are  known.  Implementation  of 
the  theory  as  a usable  algorithm  requires  estimates  of 
volcanic  load  ages  and  of  the  evolution  of  the  flexural 
response  function  with  age,  and  calculational  schemes  for 
evaluating  the  convolution  integrals  in  (5)  and  (7). 

Some  care  must  be  exercised  in  deciding  the  age  of 
a volcanic  construct  relative  to  the  age  of  the  surrounding 
seafloor.  For  islands  this  is  not  a serious  problem  as  the 
exposed  or  cored  (if  a coralline  island)  bedrocks  can  be 
dated.  For  seamounts  some  simple  rules  will  be  helpful.  One 
such  rule  is  that  an  inactive  topographic  feature  will  sink 
with  respect  to  sea  level  at  the  rate  of  its  surrounding 
abyssal  sea  floor,  a rate  which  is  associated  with  thermal 
contraction  of  the  lithosphere  and  which  is  a well  known 
function  of  ocean-floor  age  (Sclater  et  al. , 1971) . On  this 
basis,  for  instance,  the  Ninetyeast  ridge  can  be  shown  to 
have  been  generated  at  the  southeast  Indian  ridge  (Sclater 
and  Fisher,  1974).  Much  of  the  ocean  floor  topography,  in 
fact,  was  apparently  generated  at  or  near  ridge  crests  on 
very  young  oceanic  lithosphere  (McKenzie  and  Bowin,  1976) . 

The  flexural  rigidity  D and  the  effective  elastic  plate 

thickness  T are  known  to  vary  from  values  in  the  range 
27  29 

3 x 10  - 2 x 10  dyne-cm  and  3-13  km,  respectively, 

near  mid-ocean  ridges  (Watts  and  Cochran,  1977;  Detrick  and 
Watts,  1978)  to  values  in  the  range  3 x 10^  to  10^  dyne-cm 
and  15  - 3 0 km,  respectively,  for  older  lithosphere  (Walcott, 
1970;  Watts  and  Cochran,  1974;  Watts  etal.  , 1975;  1976; 

Caldwell  et  al. , 1975) . The  assumption  that  T scales  as  an 
isotherm  depth  is  consistent  with  the  data  (Caldwell  and 
Turcotte,  1978),  so  that  a standard  oceanic  plate  thermal 
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model  (Parsons  and  Sclater,  1977)  and  a few  good  measures  of 
D and  T are  sufficient  to  define  T as  a function  of  lithosphere 
age  minus  load  age. 

Evaluation  of  the  necessary  convolution  integrals  has  been 
conducted  as  follows.  Equation  (5)  is  evaluated  by  dividing 
each  area  element  over  which  mean  elevation  is  known  (e.g., 
degree  squares)  into  subelements  and  treating  the  load  on 
each  subelement  as  an  equivalent  load  on  a unit  area  in  the 
subelement  center.  Computation  time  is  saved  by  noting  that 
the  contribution  to  w(x,y)  from  the  load  q(x',y'  ) depends 
only  on  the  scalar  distance  from  (x',y')  to  (x,y) . The 
Bessel-Kelvin  functions  are  evaluated  using  a FORTRAN  sub- 
routine package  (COMBES) , written  at  the  MIT  Information 
Processing  Center,  for  calculating  Bessel  functions  of 
complex  argument  and  order.  Equation  (7)  is  evaluated 
using  the  Taylor  series  algorithm  of  Morrison  (1976)  for 
calculating  the  gravitational  potential  from  a surface 
density  distribution  specified  over  a latitude-longitude 
grid,  and  a centered  finite  difference  to  obtain  the  radial 
derivative  of  the  potential. 
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INTENDED  APPLICATIONS 


The  bathymetry-gravity  algorithm  is  currently  being 
applied  to  the  central  Pacific  basin,  a region  chosen 
because  the  gravity,  bathymetry  and  lithospheric  age  are 
all  well  known  and  because  there  are  superposed  volcanic 
loads  on  the  plate  of  various  younger  ages.  The  gravity 
for  the  region  is  taken  from  the  degree  averages  of 
Watts  and  Leeds  (1977).  Bathymetry  is  taken  from  a 
DMA  compilation  provided  by  AFGL  (T.P.  Rooney,  personal 
communication,  1978).  Plate  age  is  taken  from  standard 
magnetic  anomaly  information  (Pitman  et  ad.,  1974). 

Results  from  the  bathymetry-gravity  analysis  for  the 
central  Pacific,  and  an  assessment  of  its  general  validity, 
will  be  given  in  the  next  Scientific  Report. 
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